Imports

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math

from sklearn import linear_model
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn import tree, svm
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import r2_score
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import ExtraTreesRegressor

import torch
from torch import nn
from torch.utils.data import Dataset, DataLoader

import geopandas
import folium
from folium.plugins import HeatMap

Data

bottle1 = pd.read_csv("bottle_data_2010-14.csv")
bottle2 = pd.read_csv("bottle_Data_2015to19.csv", header=None)
cast = pd.read_csv("cast_Data_2015to19.csv")
cc = pd.read_csv('CC_Bottle_Cast.csv')
cc.shape

/Users/annaleaohalloran/opt/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3063: DtypeWarning: Columns (0,1,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,48,49,50,51,52,53,54,55,56,57,58,59,60,61) have mixed types.Specify dtype option on import or set low_memory=False.
  interactivity=interactivity, compiler=compiler, result=result)

(221588, 16)

bottle = pd.read_csv("bottle.csv", encoding='utf-8')

/Users/annaleaohalloran/opt/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3063: DtypeWarning: Columns (47,61) have mixed types.Specify dtype option on import or set low_memory=False.
  interactivity=interactivity, compiler=compiler, result=result)

print(bottle1.columns)
bottle1.head()
bottle1.info()

Index(['Cst_Cnt', 'Btl_Cnt', 'Sta_ID', 'Depth_ID', 'Depthm', 'T_degC',
       'Salnty', 'O2ml_L', 'STheta', 'O2Sat', 'Oxy_µmol/Kg', 'BtlNum',
       'RecInd', 'T_prec', 'T_qual', 'S_prec', 'S_qual', 'P_qual', 'O_qual',
       'SThtaq', 'O2Satq', 'ChlorA', 'Chlqua', 'Phaeop', 'Phaqua', 'PO4uM',
       'PO4q', 'SiO3uM', 'SiO3qu', 'NO2uM', 'NO2q', 'NO3uM', 'NO3q', 'NH3uM',
       'NH3q', 'C14As1', 'C14A1p', 'C14A1q', 'C14As2', 'C14A2p', 'C14A2q',
       'DarkAs', 'DarkAp', 'DarkAq', 'MeanAs', 'MeanAp', 'MeanAq', 'IncTim',
       'LightP', 'R_Depth', 'R_TEMP', 'R_Sal', 'R_DYNHT', 'R_Nuts',
       'R_Oxy_µmol/Kg', 'DIC1', 'DIC2', 'TA1', 'TA2', 'pH1', 'pH2',
       'DIC Quality Comment'],
      dtype='object')
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 38819 entries, 0 to 38818
Data columns (total 62 columns):
 #   Column               Non-Null Count  Dtype  
---  ------               --------------  -----  
 0   Cst_Cnt              38819 non-null  int64  
 1   Btl_Cnt              38819 non-null  int64  
 2   Sta_ID               38819 non-null  object 
 3   Depth_ID             38819 non-null  object 
 4   Depthm               38819 non-null  int64  
 5   T_degC               38819 non-null  float64
 6   Salnty               38810 non-null  float64
 7   O2ml_L               37931 non-null  float64
 8   STheta               38479 non-null  float64
 9   O2Sat                37902 non-null  float64
 10  Oxy_µmol/Kg          37898 non-null  float64
 11  BtlNum               26083 non-null  float64
 12  RecInd               38819 non-null  int64  
 13  T_prec               38819 non-null  int64  
 14  T_qual               4168 non-null   float64
 15  S_prec               38810 non-null  float64
 16  S_qual               11750 non-null  float64
 17  P_qual               4436 non-null   float64
 18  O_qual               10016 non-null  float64
 19  SThtaq               5593 non-null   float64
 20  O2Satq               7527 non-null   float64
 21  ChlorA               26230 non-null  float64
 22  Chlqua               12004 non-null  float64
 23  Phaeop               26230 non-null  float64
 24  Phaqua               12004 non-null  float64
 25  PO4uM                35415 non-null  float64
 26  PO4q                 3342 non-null   float64
 27  SiO3uM               35420 non-null  float64
 28  SiO3qu               3338 non-null   float64
 29  NO2uM                35420 non-null  float64
 30  NO2q                 3559 non-null   float64
 31  NO3uM                35418 non-null  float64
 32  NO3q                 3474 non-null   float64
 33  NH3uM                34725 non-null  float64
 34  NH3q                 4183 non-null   float64
 35  C14As1               2032 non-null   float64
 36  C14A1p               1143 non-null   float64
 37  C14A1q               34958 non-null  float64
 38  C14As2               2032 non-null   float64
 39  C14A2p               1143 non-null   float64
 40  C14A2q               34958 non-null  float64
 41  DarkAs               3019 non-null   float64
 42  DarkAp               1685 non-null   float64
 43  DarkAq               34023 non-null  float64
 44  MeanAs               3020 non-null   float64
 45  MeanAp               1685 non-null   float64
 46  MeanAq               34022 non-null  float64
 47  IncTim               2025 non-null   object 
 48  LightP               1338 non-null   float64
 49  R_Depth              38819 non-null  int64  
 50  R_TEMP               38813 non-null  float64
 51  R_Sal                38810 non-null  float64
 52  R_DYNHT              38819 non-null  float64
 53  R_Nuts               34745 non-null  float64
 54  R_Oxy_µmol/Kg        37898 non-null  float64
 55  DIC1                 1429 non-null   float64
 56  DIC2                 153 non-null    float64
 57  TA1                  1468 non-null   float64
 58  TA2                  157 non-null    float64
 59  pH1                  38 non-null     float64
 60  pH2                  4 non-null      float64
 61  DIC Quality Comment  37 non-null     object 
dtypes: float64(52), int64(6), object(4)
memory usage: 18.4+ MB

#extracting features we will use
bottle1['nitrogen'] = bottle1['NO3uM'] + bottle1['NO2uM'] + bottle1['NH3uM']
temp = bottle1[['Depthm', 'T_degC', 'SiO3uM', 'PO4uM', 'nitrogen', 'ChlorA']]
temp.dropna(inplace=True)
X = temp[['Depthm', 'T_degC', 'SiO3uM', 'PO4uM', 'nitrogen']]
y = temp['ChlorA']


X.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 24512 entries, 0 to 38818
Data columns (total 5 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   Depthm    24512 non-null  int64  
 1   T_degC    24512 non-null  float64
 2   SiO3uM    24512 non-null  float64
 3   PO4uM     24512 non-null  float64
 4   nitrogen  24512 non-null  float64
dtypes: float64(4), int64(1)
memory usage: 1.1 MB

/Users/annaleaohalloran/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  after removing the cwd from sys.path.

display(X) #inputs
display(y) #targets

	Depthm	T_degC	SiO3uM	PO4uM	nitrogen
0	0	16.360	2.80	0.32	0.000
1	2	16.360	2.80	0.32	0.000
2	5	15.920	3.00	0.34	0.000
3	10	15.620	3.30	0.37	0.120
4	11	15.600	3.30	0.37	0.120
...	...	...	...	...	...
38814	0	19.672	2.96	0.33	0.182
38815	2	19.672	2.96	0.33	0.142
38816	5	19.022	2.98	0.33	0.290
38817	10	18.722	3.06	0.34	0.199
38818	15	18.291	3.49	0.38	0.531

24512 rows × 5 columns

0        0.770
1        0.770
2        0.950
3        1.320
4        1.380
         ...  
38814    0.611
38815    0.611
38816    0.356
38817    0.486
38818    0.735
Name: ChlorA, Length: 24512, dtype: float64

y.isna().sum()

0

Exploring

Trying out different algorithms from sci-kit learn on Bottle 2010-2014

nn = MLPRegressor(random_state=1, max_iter=500).fit(X_train, y_train)
y_pred = nn.predict(X_test)
print('NN:')
print(nn.score(X_train,y_train))
print(nn.score(X_test, y_test))
print(r2_score(y_test,y_pred))

NN:
0.37286049108764974
0.34482422554367453
0.34482422554367453

rf = RandomForestRegressor().fit(X_train, y_train)
print('RF:')
print(rf.score(X_train,y_train))
print(rf.score(X_test, y_test))

RF:
0.9383785126733009
0.56755703574668

rf = RandomForestRegressor().fit(X_train, y_train)
print('RF:')
print(rf.score(X_train,y_train))
print(rf.score(X_test, y_test))

RF:
0.9395112322013726
0.5656491533978847

knn = KNeighborsRegressor().fit(X_train, y_train)
print('RF:')
print(knn.score(X_train,y_train))
print(knn.score(X_test, y_test))

RF:
0.5412531427080636
0.30504376221137286

# scaled
X_train_t= StandardScaler().fit_transform(X_train)
X_test_t = StandardScaler().fit_transform(X_test)

reg = LinearRegression().fit(X_train_t, y_train)
print('Linear regression:')
print(reg.score(X_train_t, y_train))
print(reg.score(X_test_t, y_test))

tree = tree.DecisionTreeRegressor().fit(X_train_t,y_train)
print('Decision tree:')
print(tree.score(X_train_t, y_train))
print(tree.score(X_test_t, y_test))


svm = svm.SVR().fit(X_train_t,y_train)
print('SVM:')
print(svm.score(X_train_t,y_train))
print(svm.score(X_test_t,y_test))

nn = MLPRegressor(random_state=1, max_iter=500).fit(X_train_t, y_train)
y_pred = nn.predict(X_test) # this is for r2_score
print('NN:')
print(nn.score(X_train_t,y_train))
print(nn.score(X_test_t, y_test))
print(r2_score(y_test_t,y_pred)) #same as above line

Linear regression:
0.15614814911926955
0.15701920540993064

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-13-be380a5274ad> in <module>
      8 print(reg.score(X_test_t, y_test))
      9 
---> 10 tree = tree.DecisionTreeRegressor().fit(X_train_t,y_train)
     11 print('Decision tree:')
     12 print(tree.score(X_train_t, y_train))

AttributeError: 'DecisionTreeRegressor' object has no attribute 'DecisionTreeRegressor'

Same algos but Bottle 2010-2019

#extracting features we will use
bottle['nitrogen'] = bottle['NO3uM'] + bottle['NO2uM'] + bottle['NH3uM']
temp = bottle[['Depthm', 'T_degC', 'SiO3uM', 'PO4uM', 'nitrogen', 'ChlorA']]
temp.dropna(inplace=True)
X = temp[['Depthm', 'T_degC', 'SiO3uM', 'PO4uM', 'nitrogen']]
y = temp['ChlorA']


X.info()

# goal: predict chlorA amount using inputs, and check performance by target - unscaled 
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, random_state=42)

reg = LinearRegression().fit(X_train, y_train)
print('Linear regression:')
print(reg.score(X_train, y_train))
print(reg.score(X_test, y_test))

tree = tree.DecisionTreeRegressor().fit(X_train,y_train)
print('Decision tree:')
print(tree.score(X_train, y_train))
print(tree.score(X_test, y_test))


svm = svm.SVR().fit(X_train,y_train)
print('SVM:')
print(svm.score(X_train,y_train))
print(svm.score(X_test,y_test))

nn = MLPRegressor(random_state=1, max_iter=500).fit(X_train, y_train)
y_pred = nn.predict(X_test)
print('NN:')
print(nn.score(X_train,y_train))
print(nn.score(X_test, y_test))
print(r2_score(y_test,y_pred))

rf = RandomForestRegressor().fit(X_train, y_train)
print('RF:')
print(rf.score(X_train,y_train))
print(rf.score(X_test, y_test))

knn = KNeighborsRegressor().fit(X_train, y_train)
print('KNN:')
print(knn.score(X_train,y_train))
print(knn.score(X_test, y_test))

# scaled
X_train_t= StandardScaler().fit_transform(X_train)
X_test_t = StandardScaler().fit_transform(X_test)

reg = LinearRegression().fit(X_train_t, y_train)
print('Linear regression:')
print(reg.score(X_train_t, y_train))
print(reg.score(X_test_t, y_test))

#tree = tree.DecisionTreeRegressor().fit(X_train_t,y_train)
#print('Decision tree:')
#print(tree.score(X_train_t, y_train))
#print(tree.score(X_test_t, y_test))


svm = svm.SVR().fit(X_train_t,y_train)
print('SVM:')
print(svm.score(X_train_t,y_train))
print(svm.score(X_test_t,y_test))

nn = MLPRegressor(random_state=1, max_iter=500).fit(X_train_t, y_train)
y_pred = nn.predict(X_test) # this is for r2_score
print('NN:')
print(nn.score(X_train_t,y_train))
print(nn.score(X_test_t, y_test))
print(r2_score(y_test_t,y_pred)) #same as above line

rf = RandomForestRegressor().fit(X_train_t, y_train)
print('RF:')
print(rf.score(X_train_t,y_train))
print(rf.score(X_test_t, y_test))

knn = KNeighborsRegressor().fit(X_train_t, y_train)
print('KNN:')
print(knn.score(X_train_t,y_train))
print(knn.score(X_test_t, y_test))

et = ExtraTreesRegressor().fit( X_train, y_train)
print('extra trees regressor:')
print(et.score(X_train, y_train))
print(et.score(X_test, y_test))

et = ExtraTreesRegressor().fit( X_train_t, y_train)
print('extra trees regressor, standardized:')
print(et.score(X_train_t, y_train))
print(et.score(X_test_t, y_test))

Updates so far...

• The best algorithms are: Random Forest Regressor and Extra Trees Regressor
• Unstandarzied data seems better for those algorithms(we wont use the standard scaler)
• They only reach about .61 in test scoring, which we want to be as close to 1 as possible
• Linear regression , SVM, NN and KNN give trash scores

Using CC_Bottle data

The hope here is that with this cleaner, larger data with more predictive inputs we can make better predictions and get higher scores in testing.

cc.head()

	Cast_Count	Depth	Temp	Silicate	Phosphate	ChlorA	Light_Percent	Cruise_ID	Station_ID	Date	Lat_Dec	Lon_Dec	Year	Month	Day	Nitrogen
0	20851	0	15.34	2.0	0.25	0.85	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
1	20851	1	15.34	2.0	0.25	0.85	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
2	20851	10	15.32	2.0	0.28	0.65	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
3	20851	11	15.32	2.0	0.28	0.62	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
4	20851	20	15.28	2.5	0.29	0.58	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.4

temp2 = cc[['Depth', 'Temp', 'Silicate', 'Phosphate', 'Nitrogen', 'ChlorA', 
            'Month', 'Year', 'Lat_Dec', 'Lon_Dec']]
temp2.dropna(inplace=True)
X2 = temp2[['Depth', 'Temp', 'Silicate', 'Phosphate', 'Nitrogen', 'Month', 'Year', 'Lat_Dec', 'Lon_Dec']]
y2 = temp2['ChlorA']

print(temp2.shape)

(221588, 10)

/Users/annaleaohalloran/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  This is separate from the ipykernel package so we can avoid doing imports until

X_train2, X_test2, y_train2, y_test2 = train_test_split(X2, y2, test_size=0.20, random_state=42)

rf = RandomForestRegressor(n_jobs=-1).fit(X_train2, y_train2)
print('random forest regressor training & testing scores:')
print(rf.score(X_train2,y_train2))
print(rf.score(X_test2, y_test2))

random forest regressor training & testing scores:
0.9681958216211535
0.8065068425343904

Final Model using CC_Data

• Create a null data frame with everything the same except add 2 to every val in temperature to model climate change
• Run prediction on this new dataframe using model we trained up

et = ExtraTreesRegressor(n_jobs=-1).fit( X_train2, y_train2)
print('extra trees regressor trianing & testing scores:')
print(et.score(X_train2, y_train2))
print(et.score(X_test2, y_test2))

extra trees regressor trianing & testing scores:
0.999999927608805
0.8466314375605716

prediction_data = pd.read_csv('CC_Bottle_Cast.csv')
prediction_data.head()

	Cast_Count	Depth	Temp	Silicate	Phosphate	ChlorA	Light_Percent	Cruise_ID	Station_ID	Date	Lat_Dec	Lon_Dec	Year	Month	Day	Nitrogen
0	20851	0	15.34	2.0	0.25	0.85	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
1	20851	1	15.34	2.0	0.25	0.85	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
2	20851	10	15.32	2.0	0.28	0.65	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
3	20851	11	15.32	2.0	0.28	0.62	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
4	20851	20	15.28	2.5	0.29	0.58	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.4

Prediction:

We want to predict the amount of ChlorA in the future with temperatures increased 2 degrees (simulate the effect of climate change on ChlorA amount, a proxy for phytoplankton abundance.

# add 2 degrees to each item in temp column
def add_two(x):
	return x + 2

prediction_data['Temp'] = prediction_data['Temp'].apply(add_two)

prediction_data.head()

	Cast_Count	Depth	Temp	Silicate	Phosphate	ChlorA	Light_Percent	Cruise_ID	Station_ID	Date	Lat_Dec	Lon_Dec	Year	Month	Day	Nitrogen
0	20851	0	17.34	2.0	0.25	0.85	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
1	20851	1	17.34	2.0	0.25	0.85	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
2	20851	10	17.32	2.0	0.28	0.65	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
3	20851	11	17.32	2.0	0.28	0.62	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.1
4	20851	20	17.28	2.5	0.29	0.58	NaN	1977-12-05-C-31JD	083.3 051.0	12/20/1977	33.866666	-120.141666	1977	12	20	0.4

X3 = prediction_data[['Depth', 'Temp', 'Silicate', 'Phosphate', 'Nitrogen', 'Month', 'Year', 'Lat_Dec', 'Lon_Dec']]
y3 = prediction_data['ChlorA']

predictions = et.predict(X3)
prediction_data['ChlorA'] = predictions

print(predictions)

[0.64648 0.64648 0.58838 ... 2.00075 2.16124 2.09698]

Heat Maps of Chlorophyll-A Amounts

• Existing Data: 2014-2019
• Model prediction: What will the amounts look like after we increase temperature by 2°C

# get all data per year
data2014 = cc[cc['Year']==2014]
data2015 = cc[cc['Year']==2015]
data2016 = cc[cc['Year']==2016]
data2017 = cc[cc['Year']==2017]
data2018 = cc[cc['Year']==2018]
data2019 = cc[cc['Year']==2019]
data2020 = prediction_data

d = [data2014, data2015, data2016, data2017, data2018, data2019, data2020]
titles = ['Chlorophyll-A Heat Map 2014', 'Chlorophyll-A Heat Map 2015' , 'Chlorophyll-A Heat Map 2016',
          'Chlorophyll-A Heat Map 2017', 'Chlorophyll-A Heat Map 2018', 'Chlorophyll-A Heat Map 2019' , 
          'Chlorophyll-A Heat Map After Increasing Temperatures'
         ]

hmaps = []
for i in d:
    hmap = folium.Map(location=[33.866666, -120.141666], zoom_start=7 )

    hm_wide = HeatMap( list(zip(i.Lat_Dec.values, i.Lon_Dec.values,
                            i.ChlorA.values)),
                   min_opacity=0.2,
                   radius=17, blur=15, 
                   max_zoom=1, 
                 )
    
    hmap.add_child(hm_wide)
    
    hmaps += [hmap]
    
print(len(hmaps))

7

Chlorophyll-A Heat Map 2014

display(hmaps[0])

Make this Notebook Trusted to load map: File -> Trust Notebook

Chlorophyll-A Heat Map 2015

display(hmaps[1])

Make this Notebook Trusted to load map: File -> Trust Notebook

Chlorophyll-A Heat Map 2016

display(hmaps[2])

Make this Notebook Trusted to load map: File -> Trust Notebook

Chlorophyll-A Heat Map 2017

display(hmaps[3])

Make this Notebook Trusted to load map: File -> Trust Notebook

Chlorophyll-A Heat Map 2018

display(hmaps[4])

Make this Notebook Trusted to load map: File -> Trust Notebook

Chlorophyll-A Heat Map 2019

display(hmaps[5])

Make this Notebook Trusted to load map: File -> Trust Notebook

Chlorophyll-A Heat Map Prediction

What would the map look like with an increase of 2°C?

display(hmaps[6])

Make this Notebook Trusted to load map: File -> Trust Notebook